Time dependent transport phenomena 



G. Stefanucci, S. Kurth, E.K.U. Gross 
Institut fur Theoretische Physik, Freie Universitdt Berlin, 
Armmallee 14, D-14195 Berlin, Germany, 
and European Theoretical Sepctroscopy Facility (ETSF) 

A. Rubio 

Departamento de Fisica de Materiales, Facultad de Ciencias Quimicas, UPV/EHU, 
Unidad de Materiales Centro Mixto CSIC-UPV/EHU and Donostia International Physics Center (DIPC), 
San Sebastian, Spain, and European Theoretical Sepctroscopy Facility (ETSF) 



I. INTRODUCTION 

The aim of this review is to give a pedagogical introduction to our recently proposed ab initio 
theory of quantum transport. It is not intended t o be a general overview of the field. For further 
information we refer the interested reader to Refs. I1I2C I The nomenclature quantum transport has 
been coined for the phenomenon of electron motion through constrictions of transverse dimensions 
smaller than the electron wavelength, e.g., quantum-point contacts, quantum wires, molecules, etc. 
The typical experimental setup is displayed in Fig. ^where a central region C of meso- or nano-scopic 
size is coupled to two metallic electrodes L and R which play the role of charge reservoirs. The 
whole system is initially (at time t < 0) in a well defined equilibrium configuration, described by a 
unique temperature and chemical potential (thermodynamic consistency) . The charge density of the 
electrodes is perfectly balanced and no current flows through the junction. 
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FIG. 1: Schematic sketch of the experimental setup described in the main text. A central region which also 
includes few layers of the left and right electrodes is coupled to macroscopically large metallic reservoirs. The 
system is in equilibrium for negative times. 



As originally proposed by Cini^ we may drive the system out of equilibrium by exposing the 
electrons to an external time-dependent potential which is local in time and space. For instance, we 
may switch on an electric field by putting the system between two capacitor plates far away from the 
system boundaries. The dynamical formation of dipole layers screens the potential drop along the 
electrodes and the total potential turns out to be uniform in the left and right bulks. Accordingly, the 
potential drop is entirely limited to the central region. As the system size increases, the remote parts 
are less disturbed by the junction, and the density inside the electrodes approaches the equilibrium 
bulk density. 

The Cini scheme can be combined with Time Dependent Density Functional Theory (TDDFT)£ In 
this theory, the time-dependent density of an interacting system moving in an external, time-dependent 
local potential can be calculated via a fictitious system of non-interacting electrons moving in a local, 
effective time-dependent potential. Therefore this theory is in principle well suited for the treatment 
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of nonequilibrium transport problems^ However, as far as the leads are treated as noninteracting, it is 
not obvious that in the long-time limit a steady-state current can ever develop. The reason behind the 
uncertainty is that the bias represents a large perturbation and, in the absence of dissipative effects, 
e.g., electron-electron or electron-phonon scattering, the return of time-translational invariance is not 
granted. In this review we will show that the total current tends to a steady-state value provided the 
effective potential of TDDFT is independent of time and space in the left and right bulks. Also, the 
physical mechanism leading to the dynamical formation of a steady state is clarified. 

It should be mentioned that there has been already considerable activity in the density functional 
theory (DFT) community to describe transport phenomena through systems like the one in Fig. ^ 
Most approaches are limited to the steady-state regime and are based on a self-consistency procedure 
first proposed by Lang. 7 In this steady-state approach based on DFT, exchange and correlation is ap- 
proximated by the static Kohn-Sham (KS) potential and the charge density is obtained self-consistently 
in the presence of the steady current. However, the original justification involved subtle points such 
as different Fermi levels deep inside the left and right electrodes (which is not thermodynamically con- 
sistent) and the implicit reference of non-local perturbations such as tunneling Hamiltonians within 
a DFT framework. (For a detailed discussion we refer to Ref. 0.) Furthermore, the transmission 
functions computed from static DFT have resonances at the non-interacting KS excitation energies 
which in general do not coincide with the true excitation energies. 

Our TDDFT formulation, as opposed to the static DFT formulation, is thermodynamically consis- 
tent, is not limited to the steady-state regime (we can study transients, AC responses, etc.) and has 
the extra merit of accessing the true excitation energies of interacting systems^ 

We will first use the nonequilibrium Green's function (NEGF) technique to discuss the implications 
of our approach. For those readers that are not familiar with the Keldysh formalism and with NEGF, 
in Section[n]we give an elementary introduction to the Keldysh contour, the Keldysh-Green functions 
and the Keldysh book-keeping. The aim of this Section is to derive some of the identities needed for 
the discussion (thus providing a self-contained presentation) and to establish the basic notation. In 
Section lTTTl we set up the theoretical framework by combining TDDFT and NEGF. An exact expression 
for the time-dependent total current I(t) is written in terms of Green functions projected in region 
C . It is also shown that a steady-state regime develops provided 1) the KS Hamiltonian globally 
converges to an asymptotic KS Hamiltonian when t — > oo, 2) the electrodes form a continuum of 
states (thermodynamic limit), and 3) the local density of states is a smooth function in the central 
region. It is worth noting that the steady-state current results from a pure dephasing mechanism in 
the fictitious KS system. Also, the resulting steady current only depends on the KS potential at t = oo 
and not on its history. However, the KS potential might depend on the history of the external applied 
potential and the resulting steady-state current might be history dependent. A practical scheme to 
calculate I it) is presented in Section llVl The main idea is to propagate the KS orbitals in region C 
only, without dealing with the infinite and non-periodic system^ We first show how to obtain the 
KS eigenstates tp s of the undisturbed system in Section TlV Al Then, in Section TlV Bl we describe an 
algorithm for propagating ip s under the influence of a time-dependent disturbance. The numerical 
approach of Section IIVI is completely general and can be applied to any system having the geometry 
sketched in Fig. In order to demonstrate the feasibility of the scheme we implement it for one- 
dimensional model systems in Section Here we study the dynamical current response of several 
systems perturbed by DC and AC biases. We verify that for nonintcracting electrons the steady-state 
current does not depend on the history of the applied bias. Also, we present preliminary results on 
net currents in unbiased systems as obtained by pumping mechanisms. We summarize our findings 
and draw our conclusions in Section IvTl 

II. THE KELDYSH FORMALISM 
A. The Keldysh contour 

In quantum mechanics we associate to any observable quantity O a hermitian operator O. The 
expectation (\f r |0| 1 I r ) gives the value of O when the system is described by the state \*ff). For an 
isolated system the Hamiltonian H$ does not depend on time, and the expectation value of any 
observable quantity is constant provided |\&) is an eigenstate of Hq. In this Section we discuss how to 
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describe systems which are not isolated but perturbed by external fields. Without loss of generality, 
we assume that the system is isolated for negative times t and that H(t < 0) = Hq. The evolution of 
the state is governed by the Schrodinger equation i4z\'9(t)) — H(t)\^(t)), and, correspondingly, 
the value of O evolves in time as 0(t) = (tf(i)|d|#(i)). The time-evolved state \&(t)) = S(t; 0)|*(0)), 
where the evolution operator S(t; t') can be formally written as 

In Eq. ||IJ, T is the time-ordering operator and rearranges the operators in chronological order with 
later times to the left; T is the anti-chronological time-ordering operator. The evolution operator is 
unitary and satisfies the group property S(t\ti)S(ti\t') — S(t;t') for any t\. It follows that 0(t) is 
the average on the initial state |^(0)) of the operator O in the Heisenberg representation, Ojy(i) = 
S(0;t)OS(t;0), i.e., 

0(t) = (¥(0)|S(0;t)&S(t;0)|¥(0)) = (*(0)|T e - ,i ^ d *^® 6 Te^ft d * d ®|¥(0)>. (2) 
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FIG. 2: a) The oriented contour 7 described in the main text with a forward and a backward branch between 
and t. According with the orientation the point z is later than the point z' . b) The extended oriented contour 
7 described in the main text with a forward and a backward branch between and 00. For any physical time t 
we have two points t± on 7 at the same distance from the origin, c) The generalization of the original Keldysh 
contour. A vertical track going from to — i/3 has been added and, according with the orientation chosen, any 
point lying on it is later than a point lying on the forward or backward branch. 



We can now design an oriented contour 7 with a forward branch going from t = to t and a 
backward branch coming back from t and ending in t = 0, see Fig. [21a. Denoting with z the variable 
running on 7, Eq. J5J can be formally recast as follows 

0(t) = (*(0)|T K [e-'J-r** 6 ™ 0{t)\ |*(0)). (3) 

The contour ordering operator Tk moves the operators with "later" contour variable to the left. A 
point z is later than a point z' if z' is closer to the starting point, see Fig. [2Ja. In Eq. ©, 0(t) is not 
the operator in the Heisenberg representation [the latter is denoted with Ojj{i)]. Actually, 0(t) = O 
for any t. The reason of the time argument stems from the need of specifying the position of the 
operator O in the contour ordering. 

Let us now extend the contour 7 up to infinity, as shown in Fig. [3b. For any physical 
time t there are two points z = t+ and z = t_ on 7; t— lies on the forward branch while 
lies on the backward branch and it is later than t— according with the orientation chosen. 

We have T K {e" 1 ^ dl " (i) 6(t-)} = 5(0; oo)5(oo; t)6(t)S(t; 0) = S*(0; t)OS(t; 0), and similarly 

T K {e~ l ^ diH(i) 6{t+)} = S(0;t)6(t)S(t;oo)S(oo;0) = S(0;t)OS(t;0). Thus, the expectation value 
0(t) in Eq. J2J is also given by the formula 

0{z) = (¥(0)|Zk { e - iJ "-» d ** (s) 6{z)\ |*(0)). (4) 

where 7 is the contour in Fig. [21b; 7 is called the Keldysh contourM^ In Eq. (@J the variable z can 
be either t_ or t + and 0(t_) = 0(t + ) = 0(t). 
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The Keldysh contour can be further extended to account for statistical averages^ In statistical 
physics a system is described by the density matrix p = w n \^! n )(^ n \ with w n the probability of 
finding the system in the state \H? n ) and J^n w n = 1- The states |^„) may not be orthogonal. We 
say that the system is in a pure state if p = \^)(^\. In a system described by a density matrix p(0) 
at t = 0, the time-dependent value of the observable O is a generalization of Eq. Q, i.e., O(z) = 

J2n w n(*n(0)|7k{e -1 ^-' dzH W O(z)} | (0)) . Among all possible density matrices there is one that is 
very common in physics and describes a system in thermal equilibrium: p = exp[— (3{H$ — pN)]/Z, 
with the inverse temperature (3, the chemical potential /i, the operator N corresponding to the total 
number of particles and the grand-partition function Z = Tr exp[— (3(H — pN)]. Assuming that 
Hq and N commute, the statistical average O(z) with the thermal density matrix can be written as 
O(z) = Tr [e^e-^TKje - *^ dzH(z) o( z )}]/Z. We can now extend further the Keldysh contour as 
shown in Fig. 0c and define H{z) — Hq for any z on the vertical track. With this definition H(z) is 
continuous along the entire contour since H(0) = Hq. According to the orientation displayed in the 
figure, any point lying on the vertical track is later than a point lying on the forward or backward 
branch. We use this observation to rewrite O(z) as 



Tr 


>**T K {e- i A d * A( * ) 6(z)}' 


Tr 


e /3 M JV TK | e -iX,d2ff(f)| 





0(z) = — r — 7 : : tt 1 ^, (5) 

where Xk is now the ordering operator on the extended contour. It is worth noting that the de- 
nominator in the above expression is simply Z. We have already shown that choosing z on one of 
the two horizontal branches, Eq. © yields the time-dependent statistical average of the observable 
O. On the other hand, if z lies on the vertical track O(z) = Tr [e^e^/^^Oe^^ 6 ° ]/Z = 
Tr [e~P( Ha ~^ N }(J\/Z, where the cyclic property of the trace has been used. The result is independent 
of z and coincides with the thermal average of the observable O. 

To summarize, in Eq. © the variable z lies on the contour of Fig. [2c; the r.h.s. gives the time- 
dependent statistical average of the observable quantity O when z lies on the forward or backward 
branch, and the statistical average before the system is disturbed when z lies on the vertical track. 



B. The Keldysh-Green function 

The idea presented in the previous Section can be used to define correlators of many operators on 
the extended Keldysh contour. The Keldysh-Green function G is the correlator of two field operators 
ip{r) and ip^(r) which obey the anticommutation relations {il>(r), ip'(r')} = 6(r — r'). It is defined as 



(r\G(z; 



r) = 



Tr 


e^T K 




Tr 







(6) 



where the contour variable in the field operators specifies the position in the contour ordering (there is 
no true dependence on z in ip and ifi). Here and in the following we use boldface to indicate matrices 
in one-electron labels, e.g., G is a matrix and (r|G|r') is the (r,r') matrix element of G. Due to the 
contour ordering operator Tk, the Green function G has the following structure 



G(z;z') = 6(z, z')G > {z; z') + 6(z' , z)G < {z; z'), 



(7) 



where 



[z,z 



= 1 if z is later than z 1 on the contour and zero otherwise. We say that a two-point 



function on the contour having the above structure belongs to the Keldysh space. The Green function 
G(z; z') obeys an important cyclic relation on the extended Keldysh contour. As we shall see, the 
relations below play a crucial role since they provide the boundary conditions for solving the Dyson 
equation. Choosing z — 0_ we find 



(r|G(0_;*V) 



Tr 






Tr 





(8) 
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where we have taken into account that 0_ is the earliest time and therefore V'( T *; 0_) is always moved 
to the right when acted upon by Tk. The extra minus sign in the r.h.s. comes from the contour 
ordering. More generally, rearranging the field operators ip and i/)> (later arguments to the left), we 
also have to multiply by (— l) p , where P is the parity of the permutation. Inside the trace we can 
move ip(r) to the left. Furthermore, we can exchange the position of ip(r) and e^^ N by noting that 
i/)(r)e /3AtAr = e /3/i ( Ar+1 ''0(r). Using the fact that Tk moves operators with later times to the left we 
have ij){r)Tyi{. . .} = T^{%jj(r, —iff) ■ ■ ■ }■ Therefore, we conclude that 

G(0-;z') = -c^G(-i/3;z'), G{z; 0_) = -e -ft *G(«; -i/9), (9) 

where the second of these relations can be obtained in a similar way. Eq. ||§J are the so called 
Kubo-Martin-Schwinger (KMS) boundary conditions ^iiSi 



C. The Keldysh book-keeping 

In this Section we derive some of the identities that we will use for dealing with time -dep endent 
transport phenomena. A systematic and more exhaustive discussion can be found in Ref. Ilo . 

Let k(z;z') belong to the Keldysh space: k(z;z') = 6{z, z')k > (z; z') + 9(z' , z)k < (z; z'). For any 
k(z; z') in the Keldysh space we define the greater and lesser functions on the physical time axis 

k > (t;l/) = k(t + ;t'_), k<(t,t') = ;*'+)■ (10) 

We also define the left and right functions with one argument t on the physical time axis and the other 
r on the vertical track 

fcl (t; r) = k(t ± ; r), k^r, t) = k(r; t ± ). (11) 

In the definition of k} and k^ we can arbitrarily choose t+ or t- since r is later than both of them. 
The symbols "]" and "[" have been chosen in order to help the visualization of the time arguments. 
For instance, "]" has a horizontal segment followed by a vertical one; correspondingly, k) has a first 
argument which is real (and thus lies on the horizontal axis) and a second argument which is imaginary 
(and thus lies on the vertical axis). We will also use the convention of denoting with Latin letters the 
real time and with Greek letters the imaginary time. 

It is straightforward to show that if a(z; z') and b(z; z') belong to the Keldysh space, then c(z; z') = 
dz a(z; z)b(z; z') also belongs to the Keldysh space. From the definitions Q1OI110 we find 

ft'- ft+ r—ifi 

c > (t;t') = / Az a > (t + ;z)b < {z;t'_)+ I dz a> (<+; z)b> {z;t'_) + dz a < (t + ;z)b > (z;t'_) = 

Jo_ Jt'_ Jt+ 

dta > (t;i)b < (t;t')+ / dt a > (t; *)&>(*; f) + / di a< (t; t)6>(F; t') + / df a 1 (t; f)tf(f(fy) 



Jf Jt Jo 

The second integral on the r.h.s. is an ordinary integral on the real axis of two well defined functions 
and may be rewritten as f* di a> (t;t)b> (i;f) = J°di a> (t;t)b> (t;f) + J* di a> (i)t'). Using 

this relation, Eq. (fl^|l becomes 

/•oo p — i@ 

c>{t-t')= di[a > (t;i)b A (i;t') + a R {t;i)b > (i;t')}+ df a 1 (t; f)b* (f; t'), (13) 



where we have introduced two other functions on the physical time axis 

k R (t;t') = e(t-t , )[* > (*;*') - ^{t-t 1 )], k A {t-t') = -0(1? - t)[k> {t-t') - k<(t-t')]. (14) 

The retarded function k K (t\t') vanishes for t < t', while the advanced function k A (t;t') vanishes for 
t > t'. A relation similar to Eq. I|13|l can be obtained for the lesser component c K . It is convenient 
to introduce a short hand notation for integrals along the physical time axis and for those between 
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and — 0. The symbol "•" will be used to write J °° dtf(t)g(i) as / • g, while the symbol "*" will be 
used to write J Q *^ dff(f)g(f) as / * g. Then 

c> =a> -b A + a R -b> +Q 1 c< = a< ■ b A + a R ■ b< + <J * . (15) 

Eq. (|15|l can be used to extract the retarded and advanced component of c. By definition c R (i; t') — 
9(t - f')[c > (t;i / ) - c<(t;t')} and therefore 

/'CO /"CO 

c R (i;t') =0(t-t') / dio R (t;*)[6 > (t;i / )-& < (*;*')]+^(*-0 / d * [o > (t;*)-a < (t;*)]6 A (t;* / ). (16) 



Due to the ^-function, we have f > i' for c R ^ 0. In the second term on the r.h.s. b A (i;t') contains a 
6{t' — t) and hence it must be t > i; therefore we can replace the difference in the square bracket with 
a R . Then we break the first term on the r.h.s. in two pieces by inserting ^-functions: one for t < t' 
and the other for t > t' . In compact notation we end up with 

c R = a R • 6 R , c A = a A ■ b A , (17) 

where the second relation can be proven in a similar way. It is worth noting that in the expressions 
for c R and c A no integration along the imaginary track is required. For later purposes we also define 
the Matsubara function k M (r; r') with both the arguments in the interval (0, —i/3): 

k M (T;r') = k(z = r;z' = t'). (18) 

It is straightforward to prove the following identities 

c i = a R .&i + a Ub M , J = a^ -b A + a M *tf, c m = (Z m^ 6 m (19) 

Finally, we consider the case of a Keldysh function k(z; z') multiplied on the left by a scalar 
function l(z). The function ki(z; z') — l(z)k(z;z') = 9(z, z')Z(z)fc > (z; z') + 9(z'; z)l(z)k < (z\ z') and 
hence belongs to the Keldysh space. The Keldysh components can be extracted using the definitions 
(110111114118(1 . Choosing for instance z = t+ and z' = t'_ we find kf{t; t') = /(i)fe > (t; t') and similarly 
for z — t- and z' = t' + we find kf(t;t') — l(t)fc < (t; t 1 ). More generally, the function I is simply a 
prefactor: kf = Ik*, where x is one of the Keldysh components («s, R, A, ] , |~, M). The same is true 
for k r (z; z') — k(z; z')r(z'), where r(z') is a scalar function: k* — k*r. 

III. QUANTUM TRANSPORT USING TDDFT AND NEGF 
A. Merging the Keldysh and TDDFT formalisms 

The one-particle scheme of TDDFT corresponds to a fictitious system of noninteracting electrons de- 
scribed by the Kohn-Sham (KS) Hamiltonian H s (z) = J dr&r'ip^ (r){r\H s (z)\r')ip(r'). The potential 
v s (r, t) experienced by the electrons in the free-electron Hamiltonian H s (t) is called the KS potential 
and it is given by the sum of the external potential, the Coulomb potential of the nuclei, the Hartree po- 
tential and the exchange-correlation potential v xc . The latter accounts for the complicated many-body 
effects and is obtained from an exchange-correlation action functional, u xc (r,£) = SA xc [n]/Sn(r,t) (as 
pointed out in Ref. Il7l the causality and symmetry properties require that the action functional 
A xc [n] is defined on the Keldysh contour). A xc is a functional of the density and of the initial density 
matrix. In our case, the initial density matrix is the thermal density matrix which, due to the exten- 
sion of the Hohenberg-Kohn theoremiSi to finite temperatures^ also is a functional of the density. 
We should mention here that an alternative formulation based on TDDFT has been recently proposed 
by Di Ventra and Todoro\i2£. In their approach the system is initially unbalanced and therefore the 
exchange-correlation functional depends on the initial state and not only on the density. 

The fictitious Keldysh-Green function Q(z; z') of the KS system satisfies a one-particle equation of 
motion 

■H s (z)\g(z;z') = l6(z-z'), 
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g(z;z')^-i±l-H s (z')^=15(z-z'), (20) 

with KMS boundary conditions In Eqs. I|20l) the arrow specifies where the derivative along the 
contour acts. The left and right equations of motion are equations on the extended Keldysh contour of 
Fig. 0c and S(z-z') — 4-6(z, z') — —-£p-8(z, z'). For any z ^ z' , the equations of motion are solved by 

the evolution operator on the contour S(z; z') = Tk{^ 1 ^' dz -f »( z )} 5 since i-^S(z\ z') = H s (z)S(z; z') 

and S(z; z'){— i -jp) = S(z; z')H s {z r ). Therefore, any Green function 

g(z: z') = 6(z, z r )S(z; 0_)/ >S(0_; z') + 9(z', z)S(z; 0_)/<S(0_; (21) 

with f> constrained by / > — / < = — il, is solution of Eqs. H20|) . I n order to fix the matrix / > or f K 
we impose the KMS boundary conditions. The matrix H s (z) = H s for any z on the vertical track, 
meaning that S{— «/3;0_) = e -13 ^ 3 . Eq. © then implies / < = —e^ l3( -^ a ^^f > , and taking into 
account the constraint / > — / < = il we conclude that / < = if(H s ), where f(uj) = l/fe' 3 ^ - ^- 1 + 1] 
is the Fermi distribution function. The matrix / > takes the form / > = i[f(H s ) — 1]. 

The Green function Q(z;z') for a system of non-interacting electrons is now completely fixed. Let 
us consider some Keldysh-Green functions. For z = t + and z' = we have the greater Green function 
while for z = and z' = t + we have the lesser Green function 

g > (t;t') =iS(t;0)[f(H s ) - l]S(0;t'), g < (t;t') = i S(t; 0)f(H s )S(0; t'). (22) 

Both g > and g < depend on the initial distribution function f(H s ). The diagonal matrix element of 
—ig < is nothing but the time-dependent value of the local electron density n(r, t), while ig > gives the 
local density of holes. Another way of writing —ig < is in terms of the eigenstates \ip s (0)) of H s with 
eigenvalues e s . From the time-evolved eigenstate \ip s (t)) — S(t;0)\%p s (0)) we can calculate the time- 
dependent wavefunction ip s (r,t) — (r\ip s (t)). Inserting ^ s \ip s (0)}(tp s (0)\ in the expression for g < 
we find —i(r\g < (t;t')\r') — ^ s /(£ s )'0 s (r, t)^*(r', t'), which for t — t' reduces to the time-dependent 
density matrix. Knowing the greater and lesser Green functions we can also calculate g R ' A . Taking 
into account the definitions (f 14ft we find 

g K {t- 1') = -i6{t - t')S(t; t'), g A {t; t') = i6(t' - t)S(t; t') = [g R (t' ; t)]t. (23) 

In the above expressions for g R ' A the Fermi distribution function has disappeared. The information 
carried by £/ R,A is the same contained in the one-particle evolution operator. There is no information 
on how the system is prepared (how many particles, how they are distributed, etc). We use this 
observation to rewrite g^ in terms of C/ R,A 

g$(t; f ) - g R (t; 0)£%; 0)S R (0; *')■ (24) 

Thus, g^ is completely known once we know how to propagate the one-electron orbitals in time and 
how they are populated before the system is perturbedASi For later purposes, we also observe that 
an analogous relation holds for £?1 ' ^ 

G ] (t;T)=ig R (t;0)g ] {0;r), g [ {r; t) = -i&{t- 0)a A (0; r). (25) 



B. Total current using TDDFT 



The fictitious g of the KS system will in general not give correct one-particle properties. However 
by definition g K gives the correct density n(r,t) = — i(r\g < (t;t)\r) . Also total currents are correctly 
given by TDDFT. If for instance I a is the total current from a particular region a we have 

f d f d 

I a (t)=-e J dr—n{r,t) = e J dr i— (r\g<{t; t)\r). (26) 
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where the space integral extends over the region a (e is the electron charge). We stress here that I a 
is the electronic current (the direction of the current coincides with the direction of the electrons). 

At this point, it is convenient to partition the system into three main regions: a central region C 
consisting of the junction and a few atomic layers of the left and right electrodes and two regions L, R 
which describe the left and right bulk electrodes. According to this partitioning, the KS Hamiltonian 
H s can be written as a 3 x 3 block matrix, and the left equation of motion in 1201) reads 




H LL (z) 

Hcl 




Hlc 

Hcc(z) 

Hrc 





Hcr 

H RR (Z) 



Q{z; z ) = 8(z - z')l, 



with 



G(z;z') 



Qll(z;z') Qlc(z;z') Qlr(z;z') 
Gcl{z;z') Qcc{z;z') Gcr(z;z') 
Grl(z;z) Qrc(z;z') Grr(z;z') 



(27) 



(28) 



(a similar equation is easily obtained for the right equation of motion). Choosing z on the forward 
branch of the Keldysh contour and z 1 on the backward branch of the same contour, we obtain a left 
and right equation for the lesser Green function. These equations can be used to get rid of the time 
derivative in Eq. (|26[) . We find for a = L, R 



I a (t) 



J dr (r\i-g< a (t;t)\r) 

•J dr(r\H aC g^ a (t;t)-g< c (t;t)Hc a \r) = 2eRe[TTc {Q a (t)}\, 



(29) 



where 



Q a (t) = g^ a (t;t)H a c^[G R (t,o)g < (o,o)g A (o,t) 
Gc3(t;0)Gfc(0;0)gL(0;t)H aC 



H 



Cot 



aC 



+ 



E 

f3=L,R. 



+ 



]T g cc (t;0)g< J (0;0)g A a (0;t)H aC 



~I=L,R 

E 

01=L,R 



g%(t;0)g^(0;0)g$ a (0;t)H aC 



(30) 



is a one-particle operator in the central region C and Trc denotes the trace over a complete set of 
one-particle states of C. Let us express the quantity Q a in terms of the Green function gcc projected 
in the central region. We introduce the uncontacted Green function g which obeys Eqs. (|20(l with 

H a c — Hca = 0, 



4> 



H LL {z) 







Hcc(z) 







H RR {z) 



where 



g{z;z') 



9ll{z;z') 

g cc (z;z') 




g(z;z')=6(z-z')l, 



(31) 



(32) 



9rr(z;z') 

and the same KMS boundary conditions as g. The uncontacted g allows us to convert Eqs. (|20|l into 
an integral equation which entails the KMS boundary conditions for g 

g(z-z') = g{z;z' 



= g(z;z') 



/ dzg(z]z)H oS g(z]z') 
/ dz~g(z;z)H oS g(z;z'), 

J y 



(33) 
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7 being the extended Keldysh contour of Fig. |2Jc and H g is the off-diagonal part of H s . Using the 
relations (|17|l of Section III CI we find 



/-»R,A 



»cc ' H c a g aa , 



'0a 



A 

aa ' 



(34) 



In Eq. (|30|> all matrix elements of Q < are evaluated at times (0;0). From Eq. (|15|) we see that 
c < (0; 0) = [al (0; 0), due to the theta-functions in the retarded and advanced components. There- 
fore 



g< c (0; 0) = \g\ p Hpc * Ohc] (0; 0), <?£ 7 (0; 0) 



g 1 c *ffc 7 fl.U (0;0), 



and exploiting the first two relations in Eq. Ijl9(l we also find that 

0| 7 (O;O) = S 0y g< (O; 0) + \g\ p H pa * Sec* H Cl gU (0;0). 



(35) 



(36) 



Substituting Eqs. (|34I35I36|I into Eq. and using the identities (|24I25|) for the G reen function q. 
we obtain the following expression for Q a (t) 



Qa(t) 



G R • £< 



Spa + G ■ E 



£ [ 



(*;*) 



<5/3a + G A ■ £ A 



0=L,R 

+ i 



G R (t; 0) [g 1 * • (5 Pa + G A --£i 

(f;0) 



I3=L,R 



G R -S^G f 



(*;*) 

(0;t) 



G A -£ A 



+ (G K (t; OJG^OiO)- 
where we have used the short-hand notation G = 5cc and 

£(z;z') = J2 S «' V a (z;z') = H Ca g aa (z;z')H aC 



(0;*), 



(37) 



(38) 



u=L,R 



is the so-called embedding self-energy which accounts for hopping in and out of region C. 

Having the quantity Q a (t) we can calculate the exact total current I a (t) of an interacting system of 
electrons. Eq. Ij29(l allows for studying transient effects and more generally any kind of time-dependent 
current responses. In the long time limit 



lim Q a (t) = 



G 



G • £< • G / 



(*;*) 



(39) 



provided G and S tend to zero when the separation between their time arguments increases (in this 
case, it is only the first term on the r.h.s. of Eq. (|37|l that does not vanish). This condition is 
not stringent and is fulfilled provided the electrode states form a continuum and the local density of 
states in the central region C is a smooth function. In the next Section we investigate under what 
circumstances a steady current I a develops in the long-time limit. We will also discuss the dependence 
of I a on the history of the external potential. 



C. Steady state and history dependence 

In this Section we show that a steady state develops provided 1) the KS Hamiltonian H s (t) globally 
converges to an asymptotic KS Hamiltonian when t — > oo and 2) the electrodes form a continuum 
of states (thermodynamic limit) and the local density of states is a smooth function in the central 
region. 

Let us define the asymptotic KS Hamiltonian of electrode a as H^ a — linit—joo H aa (t). The 
retarded/ advanced component of the uncontacted Green function g behaves like 

lim g* a (t,0) =ie- iH ^S, lim gt(0,t) = -iSU iH ^ (40) 

t — ► oo t— »oo 
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where S is a unitary operator and it is defined according to 



S = lim — i (41) 

t-*oo e -itl aa t 

T being the time-ordering operator. In terms of diagonalising one-body states IV'ma) °f ^™ with 
eigenvalues s^ a , the lesser component of the embedding self-energy, defined in Eq. 1)38(1 . can be 
written as 



lim £<(t;i') = lim H Cq ^ Q (i; 0) 9 < Q (0; 0)g A a (0; t')H aC 

,i — >oo t.t ' — »oo 

£ e- i [ e -t-e^a* , ]jr Oc |^~ a >^~ a |/(5Jr aa (0)5t)|^a)Wal»«o > (42) 



7 

m,m' 



where we have taken into account that g< a (0;0) = if(H aa (0)). The left and right contraction 
with a nonsingular hopping matrix H a c causes a perfect destructive interference for states with 
\ £ ma ~ £ m'a\ ~ V + an d hence the restoration of translational invariance in time 

lim V<(t;t')=iY y f ma r ma e- ie ™<* ^~ t '\ (43) 

tjt 7 — >oo * — ' 
rn 

where / mQ = (■0~ a |/(5i? QQ (O)5 t )|'0^' Q ) while r mct = i? C a I VCJ(VCJ In principle, there 
may be degeneracies which require a diagonalisation to be performed for states on the energy shell. 
The above dephasing mechanism is the key ingredient for a steady state to develop. Substituting Eq. 
(|43|l into Eq. (|39|l we obtain for the steady state current 



- 2e^/ m/3 Tr c {G R (e- /3 )r m/3 G A ( £ - /3 )Im[SA(e- /3 )]} 

m/3 

- 2e^/ ma Tr c {r ma Im[G R (e~ Q )]} (44) 

m 

with 



' elc-frgb-S^e) 



G R ' (e) = — ^gj— . (45) 



The imaginary part of G is simply given by G Im[S ]G . By definition we have 

£ R < A ( £ ) = H C RX H aC (46) 

and hence 



Im 



S R ' A (e)l = £ <^ " C)r™. (47) 



Using the above identity, the steady-state current can be rewritten in a Landauer-like^ form 

Ir = — e^2[fmL%nL — fmR%nR] = — • (48) 
in 

In the above formula T m R = J2 n ^mR an d %nL = J2 n ^mL are the TDDFT transmission coefficients 
expressed in terms of the quantities 

T"f = 2*6(6%, - e%,)Tr c {G R (e^ a )T ma G A (e^)T n0 } = T™ <*. (49) 

Despite the formal analogy with the Landauer formula, Eq. I(48|l contains an important conceptual 
difference since f ma is not simply given by the Fermi distribution function. For example, if the induced 
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change in effective potential varies widely in space deep inside the electrodes, the band structure of the 
a-electrode Hamiltonian SH aa (0)<S^ might differ from that of H^ a . However, for metallic electrodes 
with a macroscopic cross section the switching on of an electric field excites plasmon oscillations which 
dynamically screen the external disturbance. Such a metallic screening prevents any rearrangements of 
the initial equilibrium bulk-density, provided the time-dependent perturbation is slowly varying during 
a typical plasmon time-scale (which is usually less than a fs). Thus, the KS potential v s undergoes 
a uniform time-dependent shift deep inside the left and right electrodes and the KS potential-drop is 
entirely limited to the central region. Denoting with Av a (t) the difference in electrode a between the 
KS potential at time t and the KS potential at negative times, Av a (t) — v s (r £ a,t) — v s (r S a,0), 
to leading order in 1/N we then have 

H aa (t) = H aa (Q) + l a Av a (t), (50) 

meaning that H^ a = H aa (0) + l a Av^. Hence, except for corrections which are of lower order with 
respect to the system size, SH aa {0)S^ = H aa (0) and 

/ raa = /(C-AO. (51) 

The formula for the current can be further manipulated when Eq. I|51|l holds. Let us write the 
embedding self-energy as the sum of a real and imaginary part S^' A (e) = A Q (e) =F iT a (e)/2. Using 
Eq. I|47(l we can rewrite the transmission coefficients as 

%nR = Tr c ^G K (e^ lR )T mR G A (e^ lR )T L (e^ lR )^ , (52) 



T mL = Tr c {G K (e% L )T mL G A (e% L )T R (e% L )} . (53) 
Substituting these expressions in Eq. (|48(l and taking into account Eq. I|51|l we obtain 

4 S) = / ^[f(e-Av?)-f(e-Av%)}Tvc {G R (s)T L (e)G A (e)r R (s)} . (54) 

In the above equation the Green functions are calculated from Eq. I|45[) . The Hamiltonian H^ c is 
the KS Hamiltonian H s (t — ► oo) projected on region C and can be obtained by evolving the system 
for very long times. According to the Runge-Gross theorem, H^ c depends on how the system was 
prepared at t = (in our case the system is contacted and in thermal equilibrium) and on the full 
history of the time-dependent density. Therefore, the use of Eq. \5$ in the context of static DFT is 
generally not correct. Indeed, static DFT is an equilibrium theory while here we are dealing with a 
non-equilibrium process. One might argue that in the linear-response regime the static DFT approach 

is free from the above criticism. Unfortunately, this is not the case. Denoting with Sv^ the small 

(s) 

change Av^ of the effective potential in electrode a and with 5I R the corresponding current response, 
to first order in Sv^ Eq. Ij54(l yields 

<4 S) = J^T Tr c {G*ie)T , L (e)G A (e)r 0M (e)} (5v% - Svf) . (55) 

The Green functions and the T's in Eq. I|55() refer to the system in equilibrium and static DFT 
approaches can be used to evaluate the trace. However, DFT is not enough to calculate the change 
5v^ . Indeed 

6v™ - lim lhn [5v cxt (r, t) + 5V u (r, t) + 6v xc (r, t)} , (56) 

t — -too x — >±oo 

where x is the longitudinal coordinate, the plus sign applies for a — R and the minus sign for a = L. 
In the above equation v ex t is the external potential and Vh is the Hartree potential; their sum gives 
the electrostatic Coulomb potential vq, 

8v a , c =}im lim [5v CKt (r,t) + 5V H (r,t)} . (57) 
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The variation i5w xc of the exchange-correlation potential can be expressed in terms of the exchange- 
correlation kernel f xc (r,t;r' ,t') — Sv xc (r,t)/Sn(r' ,t') 



lim lim 5v xc (r,t) — lim lim / dr' / dt' f xc (r,t;r' ,t')Sn(r' ,t' 

t — >oo x — >±oo t — >oo x — >±oo / / 



(58) 



The kernel / xc depends only on the difference t — t'. We denote by /„ iXC (r', ui) the Fourier transform 
of /xc evaluated at x — ±00 for a = R, L. Then 



<5v Q ,xc = lim / ^e lut / dr'f a ^ c (r', u)8n(r', oS) 



(59) 



with 8n{r,u}) the Fourier transform of 5n(r,t). Rewriting 5v^ as 8v a .c + 8v a ^ xc and taking into 
account Eq. Ij59(l . the current response SI^ in Eq. I)55|l can also be written as 



SI 



(S) 



- / ^^Pt( £) 



(Svr.q - $v L ,c) + lim 



2n de x ' " ' t^oo J 2tt 

x / dr' (fii, KC (r',u>) - f L ^ c (r',uj))Sn(r',u)) 



(60) 

with T(e) = Tr c (Gf(e)T 0tL (s)G^(e)To lR (s)\. At zero temperature df{e)/de = <5(e-e_F), with e F 
the Fermi energy, and Eq. I|6(J|) becomes 



Wjf^ = G KS (e F ) 



(<>v R ,c - <5«l,c) + lim / ^e" 



x / dr' (f R , xc (r', w) - /i, xc (r', w)) 5n(r', to) 



(61) 



where Gks(£f) = —eT(eF)/2w is the conductance of the KS system. We conclude that also in the 
linear-response regime static DFT is not appropriate for calculating the conductance since dynamical 
exchange-correlation effects might contribute through the last term in Eq. Q61JI . Eq. (|fj If) can also 
be obtained within the framework of time-dependent current density functional theory as it has been 
shown in Ref. |2j 

We emphasize that the steady-state current in Eq. I|48() results from a pure dephasing mechanism in 
the fictitious noninteracting problem. The damping effects of scattering are described by A xc and v xc . 
Furthermore, the current depends only on the asymptotic value of the KS potential, v s (r,t — > 00). 
However, v s (r,t — > 00) might depend on the history of the external applied potential and the resulting 
steady-state current might be history dependent. In these cases the full time evolution can not be 
avoided. In the case of Time Dependent Local Density Approximation (TDLDA), the exchange- 
correlation potential w xc depends only locally on the instantaneous density and has no memory at all. 
If the density tends to a constant, so does the KS potential v s , which again implies that the density 
tends to a constant. Owing to the non-linearity of the problem there might still be more than one 
steady-state solution or none at all. We are currently investigating the possibility of having more than 
one steady state solution. 



IV. QUANTUM TRANSPORT: A PRACTICAL SCHEME BASED ON TDDFT 

The theory presented in the previous Sections allows us to calculate the time-dependent current in 
terms of the Green function Qcc = G projected in the central region. In practise, it is computationally 
very expensive to propagate G(z; z') in time (because it depends on two time variables) and also 
calculate Q a from Eq. (|37|l . Here we describe a feasible numerical scheme based on the propagation 
of KS orbitals. We remind the reader that our electrode-junction-electrode system is infinite and 
non-periodic. Since one can in practice only deal with finite systems we will propagate KS orbitals 
projected in the central region C by applying the correct boundary conditions^ 
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We specialize the discussion to nonmagnetic systems at zero temperature and we denote with 
ip s (r,0) = (r\tp s (0)) the eigenstates of H s (t < 0). The time dependent density can be computed in 
the usual way by n(r,t) = J2 occ |?/> s (r, i)| 2 , where the sum is over the occupied Kohn-Sham orbitals 
and \ip s {t)} is the solution of the KS equation of TDDFT i-^\ip s {t)) = H s (t)\i) s (t)). Using the 
continuity equation, we can write the total current I a (t) of Eq. (|26|l as 



Ia(t) 



= -e^2 f drV ■Im[ip*{r,t)'Vip s (r,t)] 

occ a 



(62) 



where n is the unit vector perpendicular to the surface element da and the surface S a is perpendicular 
to the longitudinal geometry of our system. From Eq. (|62|l we conclude that in order to calculate 
I a {t) we only need to know the time-evolved KS orbitals in region C. This is possible provided we 
know the dynamics of the remote parts of the system. As at the end of Section IIII CI we restrict 
ourselves to metallic electrodes. Then, the external potential and the disturbance introduced by the 
device region are screened deep inside the electrodes. As the system size increases, the remote parts 
are less disturbed by the junction and the density inside the electrodes approaches the equilibrium 
bulk-density. Thus, the macroscopic size of the electrodes leads to an enormous simplification since 
the initial-state self-consistency is not disturbed far away from the constriction. Partitioning the KS 
Hamiltonian as in Eq. I|27(l . the time-dependent Schrodinger equation reads 



dt 







Hll Hlc 




' \1>l) ' 






Hcl Hcc Hcr 




He) 


\1>r) 




Hue Hrr 




\*r) 



(63) 



where \ij) a ) is the projected wave-function onto the region a — L, R. C. We can solve the differential 
equation for i/jl and ipR in terms of the retarded Green function g^ a . Then, we have for a — L,R 



(64) 



Using Eq. i|64fl. the equation for ipe can be written as 



i^llMt)) = H cc (t)\i> c (t)) + I dt'S R (t,i')IV-c(0) 



HoaflL(*,0)|Va(0)), 



(65) 



a=L,R 



where S R = J2 a =L r H Ca9aa-H <*c , m accordance with Eq. 1)38(1 . Thus, for any given KS orbital we 
can evolve its projection onto the central region by solving Eq. (|65H in region C . Eq. (|65|l has also 
been derived elsewhere (for static Hamiltonians)^ To summarize, all the complexity of the infinite 
electrode-junction-electrode system has been reduced to the solution of an open quantum-mechanical 
system (the central region C) with proper time-dependent boundary conditions. 

Equation (|65|l is the central equation of our numerical approach to time-dependent transport. It is a 
reformulation of the original time-dependent Schrodinger equation (|63() of the full system in terms of an 
equation for the central (device) region only. The coupling to the leads is taken into account by the lead 
Green functions g R Q , a — L, R. Eq. Q65JI has the structure of a time-dependent Schrodinger equation 
with two extra terms. The first term describes the injection of particles induced by a non-vanishing 
projection of the initial wave-function onto the leads. The second term involves the self-energy S R 
and the wavefunction in the central region at previous times during the propagation. We will denote 
it as the memory integral. We should remark here that these memory effects are of different origin 
than those which are usually discussed in the context of TDDFT2£*2&. The latter ones arise from the 
dependence of the exchange-correlation functional on the full history of the time-dependent density. 
Most density-based functionals used at present rely on the adiabatic approximation therefore ignoring 
the functional dependence on past time-dependent densities (Ref. Wft) . 
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Equation (|65|1 is first order in time, therefore we need to specify an initial state which is to be 
propagated. We want to study the time evolution of systems perturbed out of their equilibrium 
ground state. Of course, the ground state of our noninteracting system is the Slater determinant of 
the occupied eigenstates of the full, extended Hamiltonian in equilibrium, H s (t < 0). The practical 
question then arises how one can obtain these eigenstates and how one can propagate them in time 
without having to deal explicitly with the extended Hamiltonian. Below we show how we have coped 
with these problems. 



A. Computation of KS eigenstates 

Let us consider our electrode-junction-electrode system in equilibrium (t < 0) and let ip s (r) = 
ipEj(r) be the j-th degenerate eigenstate of energy E of the KS Hamiltonian H s . The Green functions 
Q n ' A {t; t') and Q < (t; t') of the undisturbed system depend only on the difference t — t' . In absence of 
magnetic fields H s is invariant under time-reversal and the imaginary part of the Fourier transformed 
Q R is simply given by 

-hm [(r\g(E)\r')]=^S(E-E')^ E , j (r)r E ' j ( r ') • ( 66 ) 

E> j=l 

Multiplying Eq. (|66l) by ip Em (r)ipE n (i f '') and integrating over r and r' in region C we obtain 



I [ dr f drV| m (r)Im[(r|g(^)|rO]V'En(rO = E^- jB ')E' S ™i(^)^(^ , ) ) 
Jc Jc E , =1 



(67) 



where 



S m j 



(E)= f dr r Em (r)^ E j(r) = S* m (E) (68) 
Jc 



is the overlap matrix in region C between degenerate states. This matrix is Hermitian and can be 
diagonalized, i.e., 

Y^S mj (E)af\E) = X l (E)a^(E). (69) 

Next, we multiply Eq. (JBJJl by a m V (E)a i n' ) (E) and sum over m and n. The result can be written in 
terms of the new eigenf unctions aEi{r) — J2t=i a "'(^')V'Eri('") as 



- [ dr [ &r'a% l {r)Im[(r\Q{E)\r')]a EV {r') = 5 lv \ 2 l {E)Y j 5{E-E'), 
ft Jc Jc ' w 



(70) 



where we have used Eq. I|69|l and the orthonormality of the S'-matrix eigenvectors: 
J2'jZi a ^p (E)^ 1 \e) = 5ui. Equation (|7L)fl shows explicitly that the functions a_Ej(r) diagonalize 
Ini [Qcc(E)] m the central region and that the eigenvalues are positive. Since any linear combination 
of degenerate eigenstates is again an eigenstate, diagonalizing Im [Qcc{E)\ gives us one set of linearly 
independent, degenerate eigenstates of energy E. In our practical implementation described in more 
detail in Section we diagonalize 

lm[g CC (E)] (71) 



nD c (E) 



where Dc(E) = — -Tr {Im [Qcc(E)]} is the total density of states in the central region. If we use N g 
grid points to describe the central region, the diagonalization in principle gives N g eigenvectors but 
only a few have the physical meaning of extended eigenstates at this energy. It is, however, very easy 
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to identify the physical states by looking at the eigenvalues: at a given energy E only d,E eigenvalues 
are nonvanishing and they always add up to unity. The corresponding states are the physical ones. 
All the other eigenvalues are zero (or numerically close to zero) and the corresponding states have no 
physical meaning. 

The procedure described above gives the correct extended eigenstates only up to a normalization 
factor. When diagonalizing Eq. I|71[) with typical library routines one obtains eigenvectors which are 
normalized to the central region. Physically this might be incorrect. It is possible to fix the normal- 
ization by matching the wavefunction for the central region to the known form (and normalization) 
of the wavefunction in the macroscopic leads. 

It should be emphasized that the procedure described here for the extraction of eigenstates of the 
extended system from Qcc(E) only works in practice if E is in the continuous part of the spectrum 
due to the sharp peak of the delta function in the discrete part of the spectrum. Eigenstates in the 
discrete part of the spectrum can be found by considering the original Schrodinger equation for the 
full system: H s ip — Eip. Using again the block structure of the Hamiltonian this can be transformed 
into an effective Schrodinger equation for an energy- dependent Hamiltonian for the central region only: 



H cc + Hca El H HaC = E ^ cY (72) 

a=L,R a aa J 

This equation has solutions only for certain values of E which are the discrete eigenenergies of the full 
Hamiltonian H s . Since the left and right electrodes form a continuum, the dimension of the kernel 
of (E — H aa ) is zero for those energies E in the discrete part of the spectrum. We also notice that 
the second term in parenthesis in Eq. (|72|l is nothing but the real part of the retarded/advanced 
self-energy in equilibrium, see Eq. 1471) . Bound states as well as fully reflected waves will contribute 
to the density but not to the current and might play a role in the description of charge-accumulation 
in molecular transport, as, e.g., in Coulomb blockade phenomena. In our TDDFT formulation bound 
states and fully reflected waves also play an extra role, since they are needed for calculating the 
effective potential v s (which is a functional of the density) which is in turn used for extracting all 
extended states. 



B. Algorithm for the time evolution 



In order to calculate the longitudinal current in an electrode-junction-electrode system we need 
to propagate the Kohn-Sham orbitals. The main difficulty stems from the macroscopic size of the 
electrodes whose remote parts, ultimately, are taken infinitely far away from the central, explicitly 
treated, scattering region C. 

The problem can be solved by imposing transparent boundary conditions— at the electrode-junction 
interfaces. Efficient algorithms have been proposed for wave-packets initially localized in the scattering 
region and for Hamiltonians constant in time. In this Section we describe an algorithm well suited for 
delocalized initial states, as well as for localized ones, evolving with a time-dependent Hamiltonian. 

Let H s (t) be the time-dependent KS Hamiltonian. We partition H s (t) as in Section Till Bl The 
explicitly treated region C includes the first few atomic layers of the left and right electrodes. The 
boundaries of this region are chosen in such a way that the density outside C is accurately described 
by an equilibrium bulk density. It is convenient to write H aa (i), with a = L, R, as the sum of a 
term H° aa = H aa (0) which is constant in time and another term U a (t) which is explicitly time- 
dependent, H aa {t) — H° aa + U a (t). In configuration space U a {t) is diagonal at any time t since the 
KS potential is local in space. Furthermore, the diagonal elements U a (r,t) are spatially constant for 
metallic electrodes. Thus, U a (t) = U a (t)l a and C/l(£) — Un(t) is the total potential drop across the 
central region. We write H s (t) = H(t) + U(t) with 



H(t) 



H 



LL H LC 

Hcl Hcc{t) Hcr 
H RC H rr 



, and U(t) = 



U L {t)l L 

U R (t)l R 



(73) 
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In this way, the only term in H(t) that depends on t is Hccif)- For any given initial state IV'(O)) 
|?/; (0) ) we calculate \ip(t m = mAt)) = \^ m 1) by using a generalized form of the Cayley method 



l + iSH (m) 



Aijj{m) 

2 —, M {m+l) ) 
4c/ (m) 



1 _ ,'i.rr( m ) 



(74) 



with H (m) = \[H{t m+1 ) + H(t m )}, U (m) = \{U{t m+1 ) + U(t m )] and 8 = At/2. It should be noted 
that our propagator is norm conserving (unitary) and accurate to second-order in 6, as is the Cayley 
propagator. 2 ^ Denoting by \ip a ) the projected wave function onto the region a = i?, £, C, we find from 
Eq. JHl 
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iSH 
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|V4 m) ) + \S (m) ) - \M^). 



l c + iSH, 
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Here, is the effective Hamiltonian of the central region: 
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(76) 



with Hqq = ^[Hcc(t m +i) + Hcc(t m )]- The source term {S^} describes the injection of density 
into the region C, while the memory term \M^ m ^) is responsible for the hopping in and out of the 
region C . In terms of the propagator for the uncontacted and undisturbed a electrode 



9 a = 



l a + iSH° aa : 



(77) 



the source term can be written as 



2iS 



l c + iSH, 
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(78) 



with 



(to) _ 



i+^P 



and A£"'*> = n[u£>] s 



(79) 



For a wave packet initially localized in C the projection onto the left and right electrode \ipa^) vanishes 
and \S {m *>) = for any m, as it should be. The memory term is more complicated and reads 



m)\ 



l c + ioHl s ' a=L R fe=0 u a 



m—l .(m,k) 



x(irvi^)) 



(80) 



where 



Q 



(to) 



if 



Co: 



aC- 



(81) 



The quantities depend on the geometry of the system and are independent of the initial state 

^<°>. 

Below we propose a recursive scheme to calculate the Q^'s for those system geometries having 
semiperiodic electrodes along the longitudinal direction, see Fig. In this case H a a has a tridiagonal 
block form 



h a V a 
V a h a V a 
V a h a 



(82) 
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FIG. 3: Schematic sketch of an electrode-junction-electrode system with semiperiodic electrodes. 



where h a describes a convenient cell and V a is the hopping Hamiltonian between two nearest neighbor 
cells. Without loss of generality we assume that both h a and V a are square matrices of dimension 
N a x N a . Taking into account that the central region contains the first few cells of the left and right 
electrodes, the matrix Q„ has the following structure 



Q { r } = 



q { l n) 





R — 





q { ™ ] 



The qa 's are square matrices of dimension N a x N a and are given by 



la — y a 



l9 a \ 



i5H D 



V 



(83) 



(84) 



i,i 



where the subscript (1, 1) denotes the first diagonal block of the matrix in the square brackets. We 
introduce the generating matrix function 



i 



iySH 



V a , 



(85) 



i,i 



which can also be expressed in terms of continued matrix fractions 



1 



x + iySh a + y 2 5 2 V a — — 1 

x + iydh a + y 2 5 2 V a - L - V c 



V a 



-V a 
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' x + iySha + y 2 5 2 q a (x,y) 



The qa 's can be obtained from 



q ° to! 



d_ d_ 

dx dy 



q a (x,y) 



(86) 



(87) 



x=y=l 



From Eqs. JEZJ and one can build up a recursive scheme. Let us define p a 1 (x,y) — x-\- iySh a + 
y 2 S 2 q a (x 1 y) and pi" l) = + -^}p a {x,y)\ x=y=1 . Then, by definition, qi m) = V a PaV a - Using 

the identity ^[-J; + ^} m p a {x,y)p- 1 {x,y) = 0, one finds 



(i + idh a ) P ^ = (i - i8h a )pt- l) - s 2 J2(qi k) + 2 q [ t 1] + q { t 2) )p { « 



k=0 



with Pa 



(m) 
la 



for to < 0. Once q^ has been obtained by solving Eq. (|86|l with x = y = 1, we 



can calculate pa^ = [1 + i<5h a + <5 2 q'L ' ) ] 1 - Afterwards, we can use Eq. (EBl with J = Vqj4 J V Q to 



(i) 



,(1)! 



calculate and hence ; and so on and so forth. 



(i) 
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This concludes the description of our algorithm for the propagation of the time-dependent 
Schrodinger equation for extended systems. It is worth mentioning an additional complication here 
which arises for the propagation of a time-dependent Kohn-Sham equation. This complication stems 
from the fact that in order to compute |V'c™ +1 ^) a ^ time step m + 1 one needs to know the time- 
dependent KS potential at the same time step which, via the Hartree and exchange-correlation po- 
tentials, depends on the yet unknown orbitals j^^ 1 ). Of course, the solution is to use a predictor- 
corrector approach: in the first step one approximates by Hcc{tm), computes new orbitals 
l^c"^^) ano - from those obtains an improved approximation for H^q. 



V. IMPLEMENTATION DETAILS FOR ID SYSTEMS AND NUMERICAL RESULTS 

All the methodological discussion of Section Hvl is general and can be applied to all systems having 
a longitudinal geometry like the one in Fig. In this Section we show that the proposed scheme is 
feasible by testing it against one-dimensional model systems. The extension to real molecular-device 
configurations is presently under developmental! We consider systems described by the Hamiltonian 



(x\H\x') = 6(x 



(89) 



We have used a simple three-point discretization for the second derivative 



with equidistant grid points a;,, i = 1, . . . ,N g and spacing Ax. Within this approximation matrices 
of the form Hc a MH a c which are N g x N g matrices and appear, e.g., in Eq. (|38Jl or 1)81 [I. have only 
one nonvanishing matrix element. For a — L this is the (1, 1) clement, for a = R it is the (N g ,N g ) 
element. 

In order to proceed we have to specify the nature of the leads and therefore the lead Green function. 
Here we choose the simplest case of semi- infinite, uniform leads at constant potential U a Q. In this 
case, the retarded Green function g^ a in the energy domain can be given in closed form: 



[9a a ( E )]ki = — t -^=exp{i-\j2E a \x k -x l \ 



+ ~i=|p ex P Y^JzfiaQxk -x a0 \ + \xi -x a0 |)| (91) 

with E a = E — U a Q. The abscissa x a o is the position of the interface between the lead and the device 
region; in our implementation xlo is the first grid point of region C while xrq is the N g -th grid point 
of region C. According to the notation in Eq. the one-particle state of region C describing an 
electron localized in xlo is denoted by \xci) while the one-particle state of region C describing an 
electron localized in xrq is denoted by \xcN g )- The coordinate Xk = x a o ± kAx, k > 0, where the 
plus sign applies for a = R and the minus sign for a = L. 

The results of the procedure for calculating extended eigenstates as described in Section IIV Al is 
illustrated in Fig.0|for a square potential barrier with zero potential in both leads. In the left panel we 
have the square modulus of eigenstates at an energy below the barrier height while in the right panel 
eigenstates with energy higher than the barrier are shown. The states result from diagonalization of 
Eq. (|71J) . In order to obtain the normalization constant we compute the logarithmic derivative at the 
boundary of the central region numerically and match it to the analytic form in the lead to obtain 
the phase shift 5 a : 



= qcot(S a ) (92) 

X — X a 



where q = \J 2E a . Knowing the phase shift we can rescale the wavefunction such that it matches with 
the analytic form sin(g(a; — x a o) + S a ) at the interface. Of course, this form of the extended states only 
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FIG. 4: Continuum states of square potential barrier at different energies with leads at zero potential. Left 
panel: eigenstates for e = 0.45 a.u., just below the barrier height of 0.5 a.u.. Right panel: eigenstates at 
e = 0.6 a.u.. 



applies for E a > but as long as E is in the continuous part of the spectrum, it is correct at least for 
one of the leads. This is sufficient to determine the normalization. The states obtained numerically 
with this procedure coincide with the known analytical results. 

We then implemented the propagation scheme presented in the previous Section. Within our three- 
point approximation, h ai V a and q a are lxl matrices. The equation for q a ^ [see Eqs. (I8fi|) and JEZJ] 
becomes a simple quadratic equation which can be solved explicitly 



(o) _ -(1 + iSh a ) + y/(l + i5h a y + A{5V a y 
Q a — ' > 

Although the quadratic equation has two solutions, the above choice for q a ^ is dictated by the fact 
that the Taylor expansions for small S of Eqs. (|93|) and (|86J) have to coincide. Using this result we 

then solved the iterative scheme to obtain the q^ for m > 1. 

As a first check on the propagation method we propagated a Gaussian wavepacket which, at initial 
time t = 0, is completely localized in the central device region. (The source terms \S^ m ^) then vanish 
identically). As can be seen in Fig. [SJ the wavepacket correctly propagates through the boundaries 
without any spurious reflections. 

For the propagation of the extended initial states (the eigenstates of the unperturbed system) we 
also need to implement the source terms \S^). In the following we assume that the left and right 
leads are at the same potential initially so that the analytic form of the initial states is in both leads 
given by sin(q(x — x a o) + S a ) = [exp(iS a — iqx a o) exp(iqx) — ex.] /2i. Let us specialize the discussion 
to the case a — R and define the state \qu) according to (xRk\QRj = exp(iqkAx), where \xuk) is the 
one-particle state of electrode R describing an electron localized in xj~ — xro + fcAa;, k > 0. Then, the 

projection of the initial state onto lead R reads |V"jj ) = y% [ ex P(*^a)lfe) — ex P( — ^a)\ ~ Qb)]- From 
Eq. H78JI the contribution to the source term for a — R is completely known once we know how 
Hcfl[<7fl] m /(1.R + iSHjtn) acts on the state \qn}- We have 

HcR d fefe = v n\ x CN 9 )(x m \ - jf*]"' \q R ) (94) 

(1 R + ibH RR ) [1 R + ibH RR ) 

where xcn corresponds the iV a -th discretization point of region C (the last point on the right before 
electrode R starts). We rewrite the unknown quantity as follows 



lii + i5H RR ml 



(95) 

x=y=l 
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FIG. 5: Time evolution of a Gaussian wavepacket with initial width 1.0 a.u. and initial momentum 0.5 a.u. for 
various propagation times. The transparent boundary conditions allow the wavepacket to pass the propagation 
region without spurious reflections at the boundaries. 



with 



p{x,y) = (xm\- 



1 



\qr) 



xl R + iyS£ RR 

Next, we use the Dyson equation to find an explicit expression for p(x, y). We have 
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iyS 



xl R + lyoH RR x x xl R + iyoH RR 



H RR \q R ) 



It is straightforward to realize that the action of H RR on \q R ) yields 



H RR \q R ) = 2V R cos(qkx)\q R ) - V R e 



-iqA.x 



XRl) 



so that Eq. (|TJ7|) can be rewritten as 
2iySV R cos(qAir) 



xl R + iySH RR 
Projecting Eq. I|99|l on (x R q\ we find 
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(96) 
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(100) 



where q R (x, y) is the generating function defined in Eq. Ij85(l . Solving Eq. H100(l for p(x. v) we conclude 
that 



V R p(x, y) 



V R + iy8e-^ qR {x iy ) 
x + 2iySV R cos(qAx) 



(101) 
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Using the relation in Eq. (|87() for the coefficients we find 



[D{x,y)] m 

-p{x,y) 



m ! 



(1 - 2i5V R cos(qAx)) m t iS_ iqAx 



x=y=l 



(l + 2iSV R cos{qAx)) m+1 ' V B : 
^ (1 - 2i8V R cos(qAx)) m ~ : ' f (j) 
^5 (1 + 2%6 Vr cos(qAx)) v 



(102) 



One may proceed along the same lines for extracting the left component of the source term. 

To test our implementation we have propagated eigenstates of the extended system. As expected, 
these states just pick up an exponential phase factor exp(— iEt) during the propagation. 

We are now in a position to apply our algorithm to the calculation of time-dependent currents in one- 
dimensional model systems. The systems are initially in thermodynamic equilibrium. At time t = 0, 
a time-dependent perturbation is switched on. In all the examples below the current is calculated 
according to Eq. 

/ kp dfc / d 
— Im U>i(x,t)—%l> k (x,t) 

kF dfc . / , . d , , . d 



= 2 io SF^^S^ + ^diH (103) 

where the prefactor 2 comes from spin and kp = \J2sp is the Fermi wavevector of a system with 
Fermi energy ep . 



A. DC bias 



As a first example we considered a system where the electrostatic potential vanishes identically both 
in the left and right leads as well as in the central region which is explicitly propagated. Initially, all 
single particle levels are occupied up to the Fermi energy ep. At t — a constant bias is switched on 
in the leads and the time-evolution of the system is calculated. We chose the bias in the right lead as 
the negative of the bias in the left lead, Ur = — Ul- 

The numerical parameters are as follows: the Fermi energy is £p — 0.3 a.u., the bias is Ul = —U R = 
0.05, 0.15, 0.25 a.u., the central region extends from x — —6 to x — +6 a.u. with equidistant grid 
points with spacing Ax = 0.03 a.u.. The fc-integral in Eq. (|103fl is discretized with 100 fc-points which 
amounts to a propagation of 200 states. The time step for the propagation was At = 10~ 2 a.u. 

In Fig. El we have plotted the current densities at x = as a function of time for different values of 
the applied bias. As a first feature we notice that a steady state is achieved, in agreement with the 
discussion of Section IIII CI The corresponding steady-state current 1^ can be calculated from the 
Landauer formula. For the present geometry this leads to the steady current 



max((7i,C/ R ) 



%L\f( u -u L )-f(u>-u R )] 

Z7T 



y/u> - U L y/w - U R 



ul + Ur\ 2 + u L u R 



sin(^v2w) 
x/lu 



(104) 



where I is the width of the central region. From Eq. (|1()4|) with I = 12 a.u. and Ul = —Ur, the 
numerical values for the steady-state currents are 0.0316 a.u. (Ul = 0.05 a.u.), 0.0883 a.u. (Ul = 0.15 
a.u.) and 0.0828 a.u. (Ul — 0.25 a.u.). We see that our algorithm yields the same answers. Second, 
we notice that the onset of the current is delayed in relation to the switching time t — 0. This is easily 
explained by the fact that the perturbation at t — happens in the leads only, e.g., for |a;| > 6 a.u., 
while we plot the current at x = 0. In other words, we see the delay time needed for the perturbation 
to propagate from the leads to the center of our device region. We also note that the higher the bias 
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FIG. 6: Time evolution of the current for a system where initially the potential is zero in the leads and the 
propagation region. At t = 0, a constant bias with opposite sign in the left and right leads is switched on, 
U = Ul = —Ur (values in atomic units). The propagation region extends from x = —6 to x = +6 a.u.. The 
Fermi energy of the initial state is ef = 0.3 a.u.. The current in the center of the propagation region is shown. 



the more the current overshoots its steady-state value for small times after switching on the bias. 
Finally it is worth mentioning that increasing the bias not necessarily leads to a larger steady-state 
current. 
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FIG. 7: Left panel: Time evolution of the current through a double square potential barrier in response to an 
applied constant bias (given in atomic units) in the left lead. The potential is given by V{x) = 0.5 a.u. for 
5 < | a; | < 6 a.u. and zero otherwise, the propagation region extends from x = —6 to x = +6 a.u.. The Fermi 
energy of the initial state is ef = 0.3 a.u.. The current in the center of the structure is shown. Right panel: 
Time evolution of the total number of electrons in the region \x\ < 6 for the same double square potential 
barrier. 



In the second example we studied a double square potential barrier with electrostatic potential 
V{x) — 0.5 a.u. for 5 a.u. < \x\ < 6 a.u. and zero otherwise. This time we switch on a constant bias 
in the left lead only, i.e., Ur = 0. The Fermi energy for the initial state is ep = 0.3 a.u.. The central 
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region extends from x — —6 to x = +6 a.u. with a lattice spacing of Ax = 0.03 a.u.. Again, we use 
100 different A;- values to compute the current and a time step of At = 10 -2 a.u.. 

In Fig. (Left panel) we plot the current at x = as a function of time for several values of the 
applied bias U = Ul- Again, a steady state is achieved for all values of U. As discussed in Fig. El 
the transient current can exceed the steady current; the higher the applied voltage the larger is this 
excess current and the shorter is the time when it reaches its maximum. Furthermore, the oscillatory 
evolution towards the steady current solution depends on the bias. For high bias the frequency of 
the transient oscillations increases. For small bias the electrons at the bottom of the band are not 
disturbed and the transient process is exponentially short. On the other hand, for a bias close to the 
Fermi energy the transient process decays as a power law, due to the band edge singularity. As pointed 
out in Section lmCl for non-interacting electrons the steady-state current develops by means of a pure 
dephasing mechanism. In our examples the transient process occurs in a femtosecond time-scale, which 
is much shorter than the relaxation time due to electron-phonon interactions. 

In Fig. (Right panel) we plot the time evolution of the total number of electrons in the device 
region for the same values of Ul- We see that as a result of the bias a quite substantial amount of 
charge is added to the device region. This result has important implications when simulating the 
transport through an interacting system as the effective (dynamical) electronic screening is modified 
due not only to the external field but also to the accumulation of charge state in the molecular 
device. This illustrates that linear response might not be an appropriate tool to tackle the dynamical 
response and that we will need to resort to a full time-propagation approach as the one presented in 
this review. Here we emphasize that all our calculations are done without taking into account the 
electron-electron interaction. If we had done a similar calculation with the interaction incorporated in 
a time-dependent Hartree or time-dependent DFT framework we would expect the amount of excess 
charge to be reduced significantly as compared to Fig. [7| 



B. Time-dependent biases 
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In the previous Section we have shown how a steady current develops after the switching on of 
a constant bias and discussed the transient regime. Here we exploit the versatility of our proposed 
algorithm for studying different kinds of time-dependent biases. 
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FIG. 8: Same system of Fig. |7| exposed to a suddenly switched on bias at t = 0. The bias is then turned off 
at t — 75 a.u. The current is measured in the middle of the central region. 
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FIG. 9: Time evolution of the current for a square potential barrier in response to a time-dependent, harmonic 
bias in the left lead, U^it) = !7osin(wt) for different amplitudes Uo (values in a.u.) and frequency ui — 1.0 
a.u.. The potential is given by V(x) = 0.6 a.u. for \x\ < 6.0 a.u. and zero otherwise. The propagation region 
extends from x = —6 to x — +6 a.u.. The Fermi energy of the initial state is ef = 0.5 a.u.. The current at 
x — is shown. 



As a first example we consider how the current responds to a sudden switching off of the bias. For 
comparison we have considered the same double square potential barrier of Fig. subject to the 
same suddenly switched on bias, but we have turned off the bias at t = 75 a.u. The results (obtained 
with the same parameters of Fig. are displayed in Fig. |H1 We observe that the current shows a 
rather well pronounced peak shortly after switching off the perturbation. The amplitude of the peak 
is proportional to the originally applied bias. This peak always overshoots the value of the current at 
the steady state. Another interesting feature is the fact that after turning off the bias the transient 
currents show only two oscillations around the zero current limit and the transient time for switching 
off is much shorter than for switching on a high bias. 

We have also addressed the simulation of AC-transport. We computed the current for a single 
square potential barrier with V(x) — 0.6 for |x| < 6 and zero otherwise. Here we applied a time- 
dependent bias of the form Utit) = Uos'm(uit) to the left lead. The right lead remains on zero bias. 
The numerical parameters are: Fermi energy sp = 0.5 a.u., device region from x = —6 to x = +6 a.u. 
with lattice spacing Ax — 0.03 a.u.. The number of fc-points is 100 and the time step is At = 10~ 2 
a.u.. In Fig. we plot the current at x — as a function of time for different values of the parameter 
Uq = 0.1,0.2,0.3 a.u. The frequency was chosen as u> = 1.0 a.u. in both cases. Again, as for the DC- 
calculation discussed above, we get a transient that overshoots the average current flowing through 
the constriction; again, this excess current is larger the higher the applied voltage. Also, after the 
transient we obtain a current through the system with the same period as the applied bias. Note, 
however, that (especially for the large bias), the current is not a simple harmonic as the applied AC 
bias. 

Exposing the system to an AC bias also allows us to acquire information about the excitation 
energies of the molecular device. In Fig. EH (Left panel) we plot the time dependent current for a 
symmetric double square potential barrier in response to a harmonic bias in the left lead, f/z,(i) = 
Uq sin(wi), with Uq = 0.15 a.u. and u> — 0.03 a.u.. The Fermi energy of the initial state is Ef = 0.3 a.u. 
and the current at x — is shown. The central region extends from x = —6 to x = 6 a.u. with lattice 
spacing Ax = 0.03 a.u. and the potential V(x) in region C is given by V(x) = for |x|/a.u. < (6 — d) 
and V{x) = 0.5 a.u. for (6 — d) < |x|/a.u. < 6. The number of /c-points is 100 and the time step 
is At = 10~ 2 a.u.. We have studied barriers of different thickness d = 1 a.u. and d — 2 a.u.. For 
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d = 2 a.u. we observe small oscillations superimposed to the oscillations of frequency u> = 0.03 a.u. 
driven by the external AC field. Such small oscillations have frequency ~ 0.23 and can be understood 
by looking at the transmission function T(E) in the Right panel of Fig. For d = 2 a.u. both 
the second and third peaks of T(E) are in the energy window (ef — Uo,Sf + Uq) = (0.15,0.45) a.u.. 
The energy difference between these two peaks corresponds to a good extent to the frequency of the 
superimposed oscillations. On the contrary, for d = 1 a.u. only one peak of the transmission function 
T(E) is contained in the aforementioned energy window and no superimposed oscillations are clearly 
visible. This example shows the AC quantum transport can be used also for probing molecular devices. 




t/a.u. E/a.u. 



FIG. 10: Left panel: Time evolution of the current for a symmetric double square potential barrier in response 
to a time-dependent, harmonic bias in the left lead, Uh(t) = Uo sin(tjt) with Uo = 0.15 a.u. and ui = 0.03 a.u. 
for different thickness d = 1 and d — 2 a.u. of the barriers. Right Panel: Transmission function of the same 
double square potential barrier for d — 1 and d — 2 a.u. 



C. History dependence 

In Fig.^Jwe show time-dependent currents for the same double barrier as in Fig.0for two different 
ways of applying the bias in the left lead: in one case the constant bias Uq is switched on suddenly 
at t = (as in Fig. EJ, in the other case the constant Uq is achieved with a smooth switching 
U(t) — Uo sin 2 (ust) for < t < ir/(2u). As expected and in agreement with the results of Section 
IIII CI the same steady state is achieved after the initial transient time. However, the transient current 
clearly depends on how the bias is switched on. 

According to the result in Eq. (|39J) . for noninteracting electrons the independence of the history is 
not limited to steady-state regimes. The long-time behaviour of currents J(t) and I'(t) induced by 
biases U a (t) and U' a {t) does not change provided U a — U' a ^ for t — > oo. For instance, the current 
response to an AC bias has the same periodic modulation and the same phase independently of how 
the AC bias is switched on. In Fig. El we plot the time-dependent current for the same system (and 
using the same parameters) of Fig. [5J The bias remains on zero in the right lead. In the left lead 
we applied a time-dependent bias of the form Ui,{t) = Uof(t)sm(ut), with Uq = 0.2 a.u., u> = 1.0 
a.u., and we considered two different "switching on" functions f(t). The first is f(t) = 1 (as in Fig. 
EJ while the second is a ramp-like switching-on fit) = 6(T - t)t/T + 9(t - T) with T = 30 a.u.. As 
expected, and in agreement with Eq. I|39[). the current has the same behaviour in the long-time limit. 

D. Pumping current: preliminary results 

Our algorithm is also well-suited to study pumping of electrons. An electron pump is a device which 
generates a DC current between two electrodes kept at the same bias. The pumping is achieved by 
applying a periodic gate voltage depending on two or more parameters. Electron pumps have been 
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FIG. 11: Time evolution of the current for a double square potential barrier when the bias is switched on in 
two different manners: in one case, the bias Uo is suddenly switched on at t — while in the other case the 
same bias is achieved with a smooth switching U(t) — Uo sin 2 (u£) for < t < 7r/(2o»). The parameters for the 
double barrier and the other numerical parameters are the same as the ones used in Fig. |7| Uo and u given in 
atomic units. 
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FIG. 12: Time evolution of the current for a square potential barrier in response to a time-dependent, harmonic 
bias in the left lead, Uhit) — Uof(t) sin(ut) with Uo = 0.2 a.u. and frequency uo = 1.0 a.u.. The system and 
the parameters used are the same as in Fig. E] The current at x = is shown for two different "switching on" 
functions f(t). 
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realized experimentally, e.g., for an open semiconductor quantum do'b^i where pumping was achieved 
by applying two harmonic gate voltages with a phase shift. 

In the literature, different techniques have been used to discuss electron pumping theoretically. 
Brouwer 32 suggested a scattering approach to describe pumping of non-interacting electrons which 
has been used, e.g., to study pumping through a double barrier- 3 . Nonequilibrium Green's function 
techniques have been used to study pumping in tight-binding models of coupled quantum dots^i. 
Alternatively, Floquet theory which describes evolution of a quantum system under the influence of 
time-periodic fields is also well-suited to describe pumping 3 ^. 




x/a.u. 

FIG. 13: Snapshots of the density for and the travelling potential wave at various times for pumping through a 
single square barrier by a travelling wave. The barrier with height 0.5 a.u. extends throughout the propagation 
window from x — —8 to x = +8. The leads are on zero potential and the Fermi level is at 0.3 a.u.. The travelling 
potential wave is restricted to the propagation window |x| < 8 and has the form U(t) = Uo sin(gx — ut) with 
amplitude Uo = 0.35 a.u., wave number q = 1.6 a.u. and frequency u) = 0.2 a.u.. The initial density is given 
by the red line. 



As a first example of electron pumping we have calculated the time evolution of the density for 
a single square barrier exposed to a travelling potential wave U(t) = UoSXD.(qx — tot). The wave is 
spatially restricted to the explicitly treated device region which in our case also coincides with the static 
potential barrier. Some snapshots of the density and the potential wave are shown in Fig. 1131 The 
density in the device region clearly exhibits local maxima in the potential minima and is transported 
in pockets by the wave. This is also evident in Fig. 1141 where we show the time-dependent density 
as function of both position and time throughout the propagation. The density contour lines show 
transport of electrons from the left lead at x = —8 to the right lead at x = +8 a.u.. The pumping 
mechanism in this example resembles pumping of water with the Archimedean screw. 

As a second example we have calculated pumping through a double square barrier by applying two 
harmonic gate voltages with a phase difference to the barrier potentials, i.e., U(x, t) = Uosm(u)t) for 
the left barrier and U(x,t) — Uo s'm(uit + <j>) for the right barrier. Fig. 1 151 shows the DC component of 
the pump current as a function of the phase <f> which has a sinusoidal dependence for our parameter 
values. This is in agreement with similar results of Ref. |^ for small amplitudes of the AC bias which 
were obtained using Brouwer's approach. In addition, this example may be interpreted as a very 
simple model to describe the experiment of Ref. Ij^J 
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FIG. 14: Contour plot of the time-dependent density for pumping through a single square barrier by a travelling 
potential wave. The parameters are the same as for Fig. 1131 




FIG. 15: Parametric pumping through a double square barrier. The device region extends from x = —6 to 
x — +6 a.u., the static potential has the value 0.525 a.u. for 5 < \x\ < 6 a.u. and zero elsewhere in the device. 
Pumping is achieved by harmonic variation of the barriers, i.e., U(x,t) — Uo sin(uit) for the left barrier (—6 
a.u.< x < —5 a.u.) and U(x,t) = Uo sinfwi + (j>) for the right barrier (5 a.u.< x < 6 a.u.). The DC component 
of the pump current is displayed as a function of the phase <j>. The parameters are: Uo = 0.25 a.u., uj — 0.25 
a.u. and the Fermi energy is ef = 0.5 a.u.. 
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VI. CONCLUSIONS AND PERSPECTIVES 

In this review we have given a self-contained introduction to our recent approach to quantum 
transport. In essence our approach combines two well-established theories for the description of non- 
equilibrium phenomena of interacting many-electron systems. 

On the one hand there is the formalism of non-equilibrium Keldysh-Green functions. Although this 
approach in principle can be used to study interaction effects, here we only used it in the context of 
non-interacting electrons. The reason for this is that the self-energy of interacting electrons (which 
is not to be confused with the embedding self-energy) is long-range and nonlocal. In our scheme 
which partitions space in left and right leads as well as the central device region, this non-locality is 
extremely difficult to deal with in a rigorous manner. 

On the other hand, the NEGF formalism for (effectively) non-interacting electrons can easily be com- 
bined with the second approach for time-dependent many-particle systems, namely time-dependent 
density functional theory. Just as the NEGF formalism, TDDFT in principle gives the correct time- 
dependent density of the interacting system (if the exact exchange-correlation potential is used). 
Moreover, the time-dependent effective single-particle potential of TDDFT is a local and multiplica- 
tive potential which is crucial for practical use within the partitioning scheme for transport. 

In combining the NEGF and TDDFT approaches we have presented a formally rigorous approach 
towards the description of charge transport using an open-boundary scheme within TDDFT. We 
have implemented a specific time-propagation scheme that incorporates transparent boundaries at 
the device/lead interface in a natural way. In order to have a clear definition of a device region in 
Fig. ^ we assumed that an applied bias can be described by adding a spatially constant potential shift 
in the macroscopic part of the leads. This implies an effective "metallic screening" of the constriction. 
The screening limits the spatial extent of the induced density created by the bias potential or the 
external field to the central region. Our time-dependent scheme allows to extract both AC and DC 
I/V device characteristics and it is ideally suited to describe external field (photon) assisted processes. 

In order to illustrate the performance and potential of the method we have implemented it for 
one-dimensional model systems and applied it to a variety of transport situations: we have shown 
that a steady-state current is always reached upon application of a DC bias. For a harmonic AC bias, 
the resulting AC current need not be harmonic. In the case of systems at DC bias without any source 
of dissipation it is known that the steady-current is independent of the history of the process^. We 
have explicitly demonstrated this history independence for two different switching processes of the 
external bias. The history independence for non- interacting electrons not only applies for DC but 
also for AC bias which we have also demonstrated in a numerical example. Since we can compute 
current densities locally, we are not restricted to currents deep inside the leads. In one example we 
have analyzed the time evolution of the density for localized states which are only weakly coupled to 
the reservoirs. Finally, we have shown a few simple applications of our algorithm to electron pumping. 

The list of the example calculations presented here already demonstrates the versatility and flexi- 
bility of our algorithm. It includes the Landauer formalism as the long-time limit for systems under 
DC bias and allows to study transients. Moreover, it can deal with periodic time-dependent fields 
(which are usually treated with the Floquet formalism) but is applicable to nonperiodic conditions as 

Most theoretical approaches to transport adopt open boundary conditions and assume that trans- 
port is ballistic, i.e., under steady state conditions inelastic collisions are absent and dissipation occurs 
only in the idealized macroscopic reservoirs. This might be an unrealistic assumption for transport 
through single molecules, in particular when the device is not operated in the regime of small bias 
and linear response. When inelastic scattering dominates this picture is not applicable. In particular, 
experiments2£i22iSi indicate that inelastic scattering with lattice vibrations is present at sufficiently 
large bias, causing local heating of contacts and molecular devices. In addition, current-induced forces 
might lead to bond-breaking and electromigrations. 

In a joint collaboration with Verdozzi and Almbladh, one of us has included the nuclear degrees of 
freedom at a classical levels. The initial ground state (consisting of bound, resonant and scattering 
states) has been calculated self-consistently. Also, the time-propagation algorithm of Section IIV Bl 
has been generalized to evolve the system electrons+nuclei in the Ehrenfest approximation. Several 
aspects of the electron-ion interaction in quantum transport have been investigated. 

Electron correlations are also important in molecular conductors, for example, Coulomb blockade 
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effects dominate the transport in quantum dots. Short-range electron correlations seems to be relevant 
in order to get quantitative description of I/V characteristics in molecular constrictions 4 ^ 4 ^ 4 ^. In 
particular it is commonly assumed that the energy scales for electron-electron and electron-phonon 
interactions are different and could be treated separately. However, the metallic screening of the 
electrodes considerably reduces the Coulomb-charging energy (from eV to meV scale) . In this regime 
the energy scale for the two interactions merge and they need to be treated on the same footing. 
We would like to emphasize that our scheme allows for a consistent treatment of electronic and ionic 
degrees of freedom. 

It is clear that the quality of the TDDFT functionals is of crucial importance. In particular, exchange 
and correlation functionals for the non-equilibrium situation are required. Time-dependent linear 
response theory for DC-steady state has been implemented in Ref.0 within TDLDA assuming jellium- 
like electrodes (mimicked by complex absorbing/emitting potentials). It has been shown that the DC- 
conductance changes considerably from the standard Landauer value. Therefore, a systematic study of 
the TDDFT functionals themselves is needed. A step beyond standard adiabatic approximations and 
exchange-only potentials is to resort to many-body schemes based on perturbative expansions", 
iterative schemes^, or variational- functional formulations 48 . Another path is to explore exchange- 
correlation functionals that depend implicitlyS 5 ^^ or explicitly 50 .^ on the current density. 
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